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ABSTRACT 


This thesis investigates the behavior of Graphite-Epoxy 
composites subjected to fires as may occur on the decks of 
naval aircraft carriers. The analytical model consists of two 
parts; one for the determination of the temperature field 
within the composite due to a fire and the other for 
determining the stresses within the composite due to the 
temperature field. Both problems are provided one-dimensional 
finite element models. Appropriate failure criteria are 
incorporated to predict the survivability of composites in 
various fire environments. Parametric studies were performed 


with the results presented in both graphical and tabular form. 
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I. INTRODUCTION 


A. BACKGROUND 

Composite materials, due to their high strength to weight 
ratio, are becoming increasingly common as construction 
materials in today’s modern combat aircraft. Through use of 
these materials the weight to strength ratio of the airframe 
has been greatly reduced allowing for an increase in the 
respective payload per unit weight of the structure. 

Several different Navy fixed wing aircraft such as the F- 
18, F-16, F-15 and F-14 all use composite materials in their 
construction. Since these aircraft are used on board naval 
aircraft carriers, it 1S quite conceivable that these 
Structures could be exposed to flight deck fires of varying 
intensities during their service life. It 1s important for 
the Navy to assess the survivability of composite aircraft 
exposed to fires. To this end, experimental and analytical 
investigations have been undertaken. 

In 1979 testing was performed at the Naval Weapons Center 
(China Lake) to determine the relative survivability of 
Graphite-Epoxy Composite materials and 7075-T6é Aluminum Alloy 
when exposed to flight deck pool fires [Ref. 1]. These tests 
meaicatea Che susceptibility to burn through for Graphite- 


Epoxy samples was approximately the same as for the comparable 


7075-T6 Aluminum Alloy samples tested. However unlike the 
aluminum, in some cases the composite material continued to 
smolder and/or burn after application of extinguishing agents. 
In an attempt to model this smoldering/burning of Graphite- 
Epoxy Composites, Vatikiotis in 1980 [Ref. 2] developed a 
computer model which he later refined in 1982 [{Ref. 3]. This 
finite element program evaluates fiber diameters, fiber 
temperatures, air temperatures and oxygen levels throughout 
Graphite-Epoxy materials exposed to fire. The evaluations are 
with respect to time and location. 

Graphite-Epoxy material consists of high strength graphite 
fibers surrounded by an epoxy matrix. When exposed to a fire 
the epoxy matrix material burns off rather quickly at 400- 
500°F leaving only the graphite fibers supported by residual 
Combustion produces" These remaining high strength fibers 
develop thermally induced stresses due to the varying 
temperature profile across the material. If high enough 
stresses are developed, fibers may fail rendering the 
remaining material of no structural use. However, if a 
considerable amount of fibers do not fail the surviving 
material, although considerably weakened in the compression 


mode, could continue to be used in the tensile mode. 


B. OBJECTIVE OF THIS THESIS 
This thesis investigated thermally induced stress failure 


of Graphite-Epoxy material exposed to varying intensity fires. 


The analysis was carried out in two parts. In part one, the 
computer model developed by Vatikiotis in 1982 was used to 
model fires of different intensities. This program evaluated 
the temperature profiles across a porous medium consisting of 
graphite fibers as time progressed. These results were 
further utilized in the second part of the analysis. Part two 
of the analysis consisted of developing a stress model and 
finite element computer program. This stress code evaluated 
the thermally induced stresses of fibers based on the 
temperature profiles determined from the combustion code. The 
fiber stresses were calculated at each time step; failure 
Criteria applied and failed fibers removed from the stress 
analysis. Through this procedure, a quantitative evaluation 
of Graphite-Epoxy material failure with respect to time was 


obtained. 


C. THESIS OUTLINE 

Chapter II presents the physical and mathematical models. 

Chapter III develops the failure criteria used in 
determination of individual fiber failure. 

Chapter IV contains a description of the stress analysis 
program and it’s interaction with the combustion program 
developed by Vatikiotis. 

Chapter V presents case studies with results presented in 
graphical form. 


Chapter VI contains the conclusions. 


II. PROBLEM FORMULATION 


A. PHYSICAL MODEL 

The physical model assumed for this study is that of a 
Graphite-Epoxy composite plate exposed to fires of varying 
intensity as Shown in Figure 1. This plate of composite 
material could have several different nominal thicknesses 
depending on the number of laminate layers or ply assumed. 
For this study, nominal thicknesses of 0.25, 0.5, 0.75 angi 
inch are assumed, corresponding to approximately 35, 70, 105 


and 140 ply composites respectively. 


AlR FLOW 





Figure 1: The physical model under consideration. 


As the composite plate heats up it will lose its epoxy 
Matrix material between 400 and SOO°F. The matrix loss 
process is too complex to model exactly. However, for epoxy 
loss to occur, the epoxy matrix must undergo phase changes 
from solid to combustion products. These phase changes occur 
at a constant temperature throughout the material. 

This study assumes that no failure of graphite fibers 
occurs prior to the burn-off of the epoxy matrix material. The 
point at which the epoxy matrix has been reduced to residue 
will be chosen as the Starting time for all case studies. 
Because of the phase changes which occur during the matrix 
burn off process described above, the initial temperature 
profile across the remaining Graphite fiber medium is assumed 
to be a constant 400°F at the starting time. 

Air flows across the upper surface of the composite plate. 
This surface air flow produces a differential pressure across 
the thickness of the plate. At the starting time for these 
Studies, since the matrix material has been removed, this 
differential pressure will induce air flow perpendicular 


through the remaining composite plate. 


B. VATIKIOTIS COMBUSTION MODEL 

In order to obtain the temperature field required for 
stress analysis the one-dimensional finite element program 
developed by Vatikiotis was used. As seen in Figure 2 his 


combustion model assumes the composite material has been 


exposed to heat sufficient to remove the epoxy matrix leaving 
only the graphite fibers surrounded by air filled voids. The 
heat supplied by the fire is implemented through application 
of a constant heat flux boundary condition at the lower, X=0, 
Surface. Differential pressure across the plate thickness 
forces air flow through the material. Prior to fibem 
combustion the incoming air is considered to be of standard 
composition containing no fire combustion products. After 
combustion begins carbon dioxide and carbon monoxide also flow 


through the material. 


Graphite Fibers 


AIR FLOW 


Pressure 
Differential 


CONSTANT HEAT FLUX 





Figure 2: The Vatikiotis combustion model. 


In his presentation Vatikiotis divides the material across 


its thickness into 30 equal thickness elements having a total 


Seo) nodal points. He further develops a system of four 
transient one-dimensional equations consisting of an energy 
balance on the graphite material, an energy balance on the 
iacernal “air, conservation of species on oxygen and a 
combination of Darcy’s law and continuity equation for 
internal air flow velocities. 

The two energy balances considered a differential volume 


of material with the following convention: 


meat THEO Heat Dees Ge ois Increased 
dv Generation dV Internal Energy 


The energy balance for the graphite fibers accounted for 
heat transfer into and out of the differential volume by three 
modes; conduction, convection and radiation, using first term 
Taylor series approximations. The heat generation term, 
resulting from burning of the graphite fibers, was modeled 
using an expression of Arrhenius type. 

The energy balance for the internal air considers only 
conduction and convection at the differential volume surfaces 
along with energy transport due to the flow of air through the 
volume. Also this equation does not contain a heat generation 
werm. 


The oxygen mass balance started from the basic form of: 


ea gen ineo Oxygen OUL OF | Oxygen fi Oxygen - 
dv dV Consumption Accumulation 


The oxygen into and out of the differential volume 


considers both molecular diffusion and convective transport 


across the differential volume surfaces, again using a first 
term Taylor series approximations. The oxygen consumption was 
based on the graphite fibers (carbon) burning with the oxygen 
present, in a appropriately specified stoichiometric ratio. 
This reaction is assumed to produce carbon dioxide and carbon 
monoxide only. A relationship for the consumption of oxygen 
was developed through combination of the stoichiometric ratio 
with the Arrhenius reaction rate equation. 

The final equation of this system is the equation for 
internal air flow velocities. This equation was developed 
using a combination of Darcy’s law for flow in a porus medium 
and the continuity equation. 

The above four equations forma set of coupled, nonlinear, 
transient, partial differential equations. These equations 
are solved for graphite temperature, alr temperature, oxygen 
concentration, and internal air flow velocities with respect 
to time and position. Incorporated inthis mathematical model 
are the variations of temperature dependent properties namely, 
graphite thermal conductivity, fluid (air) viscosity, and air 
density, along with parameters derived from them. These 
properties and parameters are recalculated at each time step 
Gf -the “transient Solutsaca. The Darcy’s law-continuity 
equation was solved using a shooting method and the remaining 
three equations were solved using a Galerkin formulation of 
the finite element method. For completeness the above 


equations in their final form along with their associated 


boundary conditions are contained in Appendix A. For more in 
depth review the reader is referred to Vatikiotis’ original 
works. 

The main observations of Vatikiotis’ earlier work are very 
obvious when his computer code is used. The differential 
pressure across the porous medium causes a flow of air to flow 
through medium. This induced air flow provides for both 
convective heat transfer of the fibers and a ready supply of 
oxygen for the combustion process. These two effects of the 
air flow are conflicting effects. First, the supply of oxygen 
increases the rate of combustion which increases the fiber 
temperature. Second, the convective heat transfer away from 
the hot fibers reduces the fiber temperature. A system 
comprised of a plate and a fire is defined by (a) the fire 
intensity and its time duration and (b) the plate subjected to 
the fire. A plate with known initial temperatures and oxygen 
concentrations, thickness, porosity, and pressure differential 
is subjected to a fire of specified intensity for a given time 
period. When the fire is removed the plate may either 
continue to burn (plate combustion) or the plate may return to 
ambient conditions (plate extinction). Which process occurs 
depends on the state of the plate when the fire has been 
removed. If the temperature and air velocity fields exceed a 
Critical state then the convective heat transfer from the 
fibers to the air is overtaken by the heat generated by the 


merDer Combustion and the material will continue to burn. If 


the critical state was not achieved, the convective sean 


transfer will predominate resulting in cooling and extinction. 


C. STRESS ANALYSIS MODEL 

The combustion model described previously determines the 
temperature distribution across the thickness of the medium. 
The stress analysis model uses the temperature profiles from 
the combustion code to determine the thermally induced 
Stresses in the fiber composite plate. 

As seen in Figure 3 the stress model was developed 
assuming the graphite fibers were cylindrical rods running 
parallel to the plate surfaces. These fiber rods are modeled 
as rigidly fixed at one end and on the opposite end are 
modeled to move in unison with the other fibers. The 
expansion or contraction of fibers due to varying temperatures 
will cause forces to develop in the other fibers (See Figure 
4). These resulting forces could cause excesSive stresses to 
develop possibly causing failure. The final solution to this 
Statically indeterminate problem required use of finite 


element techniques. 


10 





Figure 3: The stress analysis model. 
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Figure 4: The effect of heating graphite rods. 
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POSITION 


As seen in Figure 5 each element (fibers per ply) was 
modeled as a linear elastic spring with forces applied at each 
end (nodal point). The finite element "direct approach" [Ref. 
7} was applied to each individual element to obtain the 
element stiffness matrix and force vector as follows. All 
forces, displacements, stresses, strains and element 


Characteristics are in the axial direction only. 





Figure 5: The linear spring element of the finite element 
direct approach. 


First referring to Figure 5, using summation of forces: 


)) F=F,\+F,=0 635) 


Ie pe ( Zoi) 
where: 


ie eOr ec applied (ae point 1 


rT] 
N 
II 


force applied at point 2 


3 


Assuming a linear spring model yields: 
Pk Cae) (2. 


where: 
u, = displacement omgeorne 1 
uy = displacement of point 2 
k = stiffness coefficient. 
From [Ref. 4] the stiffness coefficient for a cylindrical 


rod with axial displacement is: 


a (229) 
where: 
A = cross sectional area of the bar (fibers per ply) 
E = elastic modulus of fiber (axial) 
i, = lengtaset £1bem 


Combining equations (2.1a) and (2.2) yields the element 


| : 4 4 7 
n (U2), Ales 


for fibers, where n is the ply number? The forces | oo 


matrix equations of: 


K =< 
a ee 





are equal in magnitude and come from the thermal expansion of 


the material as follows. 
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The thermal strain of the material is defined as: 
Steeda ncaa Tac) CZeo) 


where: 


€-, = thermal strain 


£2 
II 


thermal coefficient of expansion 


a the temperature of the material 

T, = the temperature at which the strain is zero. 
Multiplication of equation (2.5) by Youngs modulus E yields 
the stress o, 


o=Fe=Ea (T-T,) (2G) 


which results if free expansion or contraction is prevented. 
Multiplying equation (2.6) by the cross sectional area of the 


element yields the element forces: 
ne al Fe Yl 0 ee Cu) 


at the ends of the elements where n is the ply number. 
Dropping the subscripts 1 and 2 on the forces and combining 


equation (2.7) with (2.4) and (2.3) gives the following matrix 


| ze, (2.8) 
n (42), PY 


k = AE/L based on fiber properties and dimensions 


equations for each element: 


oi 
-K_ kK 





where: 


moe — Gicplacements Of: poImes Bee 2 for fiber n 


B= SAE (T-T.) 


1 


Using the numbering system of Figure 6, the global 


stiffness matrix and global force vector are assembled as: 


(258) 
k, -k, OC g:0 OneeOeiOmniCl 0 meme F, 
-k, ky 0  Qeee0 0 | OpecQusGuiolte) -ulmeom F, 
0 0 & -k, 0 OG 6 90 00 Chememlnna F, 
0 0 -k, ky 0 “OOo cmomene: F, 
0 0 GO Oho uke ij 9Oi - CRNCRNONRCICS )iy Aiuuny mes 
0 0 0 0 -kj kj 0 0 ooeerlly| le 
0 60 0 0 O 0 Kk, =k. 0 cmememma F, 
6 0 0 0 © OO -K, ck, © Oates F, 
Popo Popo 





where: 
k, = a component in the global stiffness matrix 
U, = a component in the global displacement vector 
EF, = a component in tne gleball Fores) Veccorm 


The boundary conditions imposed on the above system of 


equations are as follows. For the rigidly fixed nodal points, 
U, =U,=U,=U,=,,,.=U,,.,=0 (2) ren 
and for the nodal points that move in unison, 


U,=U,=0,=0,=....=9; 


C2 leds 


n 

The system of equations of equation (2.9) with the 
boundary conditions of (2.10) and (2.11) applied were solved 
for nodal point displacements using IMSL math library 


Subroutine DLSARG. 


16 


Elements 





Figure 6: The numbering sequence for the finite element 
SOlwE LON . 


The stresses of the fibers per ply were calculated using 
the displacements obtained from the above system of equations 


as follows (see figure 4): 


Cece, (Zee) 


€n~€a” Eth 


m 
il 


a = strain from actual motion of nodal points=(u2-u,)/L 


m 
if 


, = Strain causing fiber stress 


€-, = fiber strain for free expansion=a(T-T,). 


sey 


Substitution of appropriate expressions for ¢, and €,, gives 


an expression for mechanical strain of: 


cnn €e Ene a (77 ; (2.13) 


Finally multiplying equation (2.13) by Young’s modulus yreldee 
(Up -u,) 
go = Fe, = Ey) -@ (T-T,) (2 , 545) 


where go is the stress developed in the fiber. 

Using the preceding equations a computer stress program 
was developed utilizing the Vatikiotis combustion code for 
fiber temperature profiles. The stress program determined 
fiber stresses throughout the composite material as a function 
Of time cand location: The stress program also determined 
failed fibers based on failure criteria of Chapter III. 


Specifics of the stress program are contained in Chapter IV. 
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III. FAILURE CRITERIA 


1 INTRODUCTION 

The stress analysis code determines fiber stresses 
throughout the composite plate and then determines if failure 
of any of the fibers has occurred. In this chapter possible 
fiber failure criteria are considered. The four failure modes 
considered were tensile strength failure, compressive strength 
mbure, Liber buckling, and fiber consumption. Of the four 
modes considered only three were incorporated in the stress 


analysis code. 


B. TENSILE FAILURE 

The first criteria considered was failure due to exceeding 
the fibers tensile strength of 300 ksi. When subjected to 
heat graphite fibers tend to contract, which is reflected in 
their thermal expansion coefficient of 20), Sspclie, 2//Shoe 
Initially, the fiber closest to the heat source will increase 
in temperature then over time the remaining fibers will heat. 


The plate fibers will continue to rise in temperature as long 


as the heat source is applied. Since thermal strain is 
proportional to temperature, the magnitude of fiber 
contraction increases with fiber temperature. Due to the 


fibers closest to the fire being hotter than those further 


from the fire, fibers closest to the fire would, if free to do 


eS 


so, contract more than those farther away. This variation in 
tendency to contract causes induced stresses since all fibers 
are restrained to contract uniformly as seen in Figure 4. 

In order to estimate the temperature required to reach 
tensile strength, a limiting case of Equation (2.14) was used. 
Since the fiber closest to the heat source is the first to 
heat, its temperature begins to rise before any significant 
change in the other fiber temperatures. As a worst case, this 
fiber could be heated to a very high temperature while all the 
remaining fibers remained very close to their initial 
temperature. If the number of fibers were relatively high, 
this situation could be modeled as one heated fiber with fixed 
ends (i.e. not allowed to move axially). Using this model 


equation (2.14) becomes: 
SS Ee ee) (3) 


or, by rearranging, the minimum temperature required for this 


fiber to reach the tensile strength is: 


S) 
Tee Seat be (3-2 





where: 
S+ = tensile strength of the fiber (300 ksi) 
Tai, = the minimum temperature for tensile failure 
E = elastic modulus (31.0X10® psi) 

T, = reference temperature (assumed to be 60°F) 


@ = thermal expansion coefficient (-0055X106e7 “ae 
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Solving the above equation for the minimum fiber 
temperature gives a value of approximately 17,600°F. This 
temperature 1S unachievable since graphites melting point is 
approximately 6000°F and, moreover, the combustion code is 
unable track temperatures over approximately 2600°F. 

Due to the above calculation, the probability of any 
fibers failing in the tensile mode is highly unlikely. 
However, a quick check for tensile failure was incorporated 
into the computer program as verification of the above 
calculations. Tensile failure was not observed during any of 


the case studies. 


C. COMPRESSIVE STRENGTH FAILURE 


The next criteria considered was purely axial compression 


of the fibers. For the graphite fibers under consideration 
their compressive strength is 260 ksi. However, this high 
value is quite deceiving. These fibers have a very small 


cross section and relatively long length, causing them to be 
Susceptible to buckling failure when in compression. Stresses 
at buckling are much lower than the compressive fracture 
Strength. For this reason compressive strength of the fibers 


weasenot used as a failure criteria. 


D. FIBER BUCKLING 
The third failure mode considered was fiber buckling. The 


graphite fibers have a very high slenderness ratio due to 


Zak 


their small cross sectional diameter and relatively lene 


length. These fibers could be susceptible to Euler column 
bulieklingr Evaluation of this mode of failure requires a 
"typical" length for the unsupported fibers. This fabes 


length was assumed to be 10 inches, corresponding to a 
reasonable sized hole burned through the material. These 
fibers were further assumed to have fixed end conditions. 
The development of Euler Buckling formula is available in 
several texts and will not be redeveloped here. From [Ref. 5] 
the formula for Euler column buckling assuming fixed end 
conditions with constant circular cross sectional area is: 


Sone 
lie 


cr 


(3.29 





where: 


JO.» = Critical stress tow column Duce 


E = elastic modulus of the fibers 


D the Diameter of the fibers 


L the length of the unsupported fibers. 


Using the above formula, the calculated value of (¢peaee 
approximately 19 psi. This value is very small and quite 
conservative since the model assumes no lateral support. 

As the epoxy matrix burns, some of the combustion products 
will stay behind forming a new matrix material. This new 
matrix structure consisting of residual combustion products, 
although weaker than the original epoxy matrix, renders 


lateral support to the fibers. This additional support 


PRP 


Significantly increases the critical buckling strength of the 
fibers. 

Roccimeret. 6). 1m his work on Gritical buckling stress of 
composite materials, develops two analytical models for 
buckling of composites, namely; the extensional mode and the 
shear mode. Of these two modes, the mode most likely for 
fiber volume fractions of approximately 20% is the shear mode. 


The shear mode formula is: 


G 
¢-..= ——— 2 (3.4) 
as Ve(1-V,) 
where: 


Ce ee oct (emle lea wOUCKEIAgGuseress 


fiber volume fraction 


Ve 


Gre (shear Modul usmot Ehe matrix material . 

A search of literature provided little information on the 
properties of epoxy combustion products for use as a matrix 
material. However, after some consideration, this problem was 
overcome by assuming that the effect of the combustion 
products, aS a matrix material, could be modeled as a small 
fraction (F) of the original graphite-epoxy composite buckling 
critical stress given by equation (3.4). This fraction (F) is 
associated with the amount of lateral Support the graphite 
fibers receive from the combustion products. Using this 
assumption the following equation was used for fiber buckling 


@riteria, 


PAE 


where oe is the reduced buckling critical stress and shame 
the combustion product fraction, that is, the fraction Ofmen= 
original stress due to combustion product lateral Sstppomee 
Failure of the graphite fibers due to exceeding the reduced 
buckling stress as calculated above is used by the stress 
analysis program. Different values of the combustion product 
fraction (F) are assumed to assess the effect of this fraction 


on failure. 


E. FIBER CONSUMPTION 

Fiber consumption was considered as the last failure mode. 
As the plate burns, under some conditions, the supply of 
oxygen in the interior of the plate becomes depleted. At that 
time the burning process moves to the X=0 surface where the 
alr is being drawn into the plate. From then on the plate 
continues to burn at the surface of the plate slowly 
decreasing in thickness as the temperature profile across the 
plate tends to a constant value. This mode of burning is 
surface recession. 

The relatively constant temperature profile produced 
across the plate in the surface recession mode, could reduce 


the stresses seen by the fibers. If these stresses are small 
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in magnitude, the possibility of fiber failure by consumption 
is possible. Surface recession fiber consumption was 
incorporated into the stress analysis code to eliminate burnt 


away fibers from the analysis. 


Z> 


IV. COMPUTER PROGRAMS 


A. MODIFICATIONS TO THE COMBUSTION CODE 

Only superficial modifications were made to Vatikiotis’ 
combustion code for this study. First, the computer code was 
brought up to Fortran 77 language by changing a few read and 
write format statements. Second, a single statement reading 
numerous input variables from the input file was broken up 
into individual read statements in order to facilitate an easy 
to read comment type input file format. This input file 
containing a short description of each input variable is 
contained in Appendix D. Third, a few output statements were 
added to facilitate building output data files for use in the 
stress analysis. The fourth and final change to this program 
was to increase the number of nodal points so that the present 
program could handle up to 101 nodal points vice the older 31 
nodal pointes. Increases in computer speed along with 
Significant memory resources over the past few years 
facilitated this change. The amount of accuracy was 
Significantly improved by using a larger number of nodal 
points, as seen in both Figure 8 and Table 1. Both display 
the same typical run of the combustion code varying only the 
number of nodal points used. Temperatures at the locations of 


X=0" (the fire surface), X=0.25" (half the medium thickness), 
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and X=0.5" (the surface opposite the fire) are compared. As 
seen, previous models of the program were not grid independent 


for the X=0" location. 


NODAL POINT EFFECT ON TEMPERATURE 
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Figure 7: A comparison of the combustion code number of 
nodal points and calculated temperature. 
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TABLE 1: INCREMENTAL TEMPERATURE CHANGE VERSUS NODAL POINTS 
USEE? 


TIME OF RUN: 60 seconds 

HEAT FLUX APPLIED: 10,000 BTU/ft? hr 
DIFFERENTIAL PRESSURE ACROSS THE MEDIUM: 0.8 psi 
MEDIUM THICKNESS: 0.5 inches 


GRAPHITE TEMPERATURE °F AT LOCATION 


POINTS 
USED X=0 inch X¥=0.25 inch X¥=0.5 inch 


822.46 607.35 | age. day 
960.57 625.83 486.23 


INCREMEN 4.4% 3.0% 0. 05% 
cutee 


1066.23 peo 493.31 
1068.11 493.14 


INCREMEN 0.2% 0.02% 0.03% 
ae 


100 Picco cop! | 672.21 | Ze aor: . 42 


101 jaO' C3 672.43 494.54 


INCREMENTAL 0.07% | ere || 0.02% 
CHANGE 
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B. STRESS ANALYSIS PROGRAM 

The stress program, as shown in Figure 8, obtains 
temperature information from an output data file produced by 
Vatikiotis’ combustion code. Location data from the 
combustion code was presented in normalized (X/L) form since 
the equations Vatikiotis used were non-dimensionalized. [In 
the stress code these locations were converted into absolute 
iecation points ({X). Until burn off of the fibers starts 
(l.e€. prior to surface recession) the non dimensionalized 
locations (X/L) remain fixed. However, when surface recession 
Starts, combustion moves to the X=0 surface and the thickness 
of the material decreases with time. As the thickness of the 
material changes the stress code continually verifies stress 
code fiber locations removing burned away fibers from the 
stress analysis. 

Vatikiotis in his combustion code formulation used a 
method of "Smearing" to account for both the fiber and air 
filled void as shown in Figure 9. Denoting the porosity of 
the medium as p, we have dV,=pdV and daV,=(1-p) dv where dV 1s 
mie fetal differential volume of air and graphite fiber, dv. 
is the differential volume of air and dV, is the differential 
volume of graphite. The effect of this smearing process is to 
replace individual fibers and voids, which are extremely large 
in number, by a homogenous material. This process is utilized 


in the basic differential equations presented in Appendix A. 
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Figure 8: Flow chart of the Stress Analysis Code. 
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"“SMEARING” PROCESS USED 
IN THE COMBUSTION CODE DEVELOPMENT 


GRAPHITE 


p dV (1-p) dV 


Where: p ia the porosity (i.e. fraction of air) 





Figure 9: The process of "Smearing" used in the combustion 
code. 


The stress code accommodated the idea of a homogeneous 
material and utilized temperature output by the following 
procedure (see Figure 10). First the graphite temperatures 
along with nodal point locations are determined from the 
combustion code. These combustion code locations are compared 
to the center location of each ply. When the ply center 
location is bracketed by the combustion code locations the 
stress code performs quadratic interpolation to determine the 
graphite fiber temperature for that ply. By this method, 
accurate temperature profiles can be obtained from the 
combustion code using 101 nodal points while at the same time 


the number fibers used in the stress code may be varied. 
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Interpolation method from 


temperatures to stress code temperatures. 


Figure 10: 
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Elemental matrix equations are determined using 
Equation (2.8). These elemental matrices are assembled into 


a global set of equations of the general form: 


fe) | Ag ee I ee hee U; FP, 
k Ky: ee oe Lae Gas Mee . | ae ine 
He : 
m | 

Kn Kas ~ Kyo Ke a B 


where: 

K = the global stiffness matrix 

U = the global displacement vector 

F = the global force vector 

m = the total number of nodal points (2n) 

j and k = are arbitrary integers where i<j and k«m. 
The boundary conditions imposed on this matrix equation (4.1), 
eee, two different forms; Type I and Type II. Their 
application is performed as follows. The Type I boundary 


Soenatti0oOns consist of those that fit the form: 


Where: 
U; = displacement at node j, 
b 


= a known constant. 


oe 


This type of boundary condition is directly comparable to the 
fixed end of the fibers (i.e. the locations of all the odd 
numbered nodal points) where b=0. Applying a single boundary 
Cone tt1o0n, U;=b to the equation set (4.1) above gives the 


following modified matrix equation: 


sl J kK m 
| 
i Kee ON Re oe on 1d eer Gi 
I +| 0 1 0 0 U;| = b (4.3) 
pig OP Re ||) Oe FE 25:5 
m ee ee OV Oo U., Fo peas 


The stress code was written to apply each fixed end boundary 
condition using the same method as (4.3) above. 


The Type II boundary conditions fit the form of: 
Op eC OS AS (4.4) 


where: 
U; and U, = are two degrees jorerrecedem 
a = known constant 


Cc 


known constant. 
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Application of the general boundary condition of Equation(4.4) 


to the matrix Equation (4.1) above yields: 


(4.5) 
a a k m 
a Ky, 0 (K,,- aX, ;) Kim U; Py -CS yy 
Sm Sp1-2Ks2) = +2 + (Ky -2akjgta2K,,+a7)- (K,,-aXj_)| | Uy By @F; CK, ack,,;*ac 
| Ke, ~ Oo» (Kg 7 aK ns) z Kn Un| | Pry CK nj 








The detail for developing this modified matrix equation is 
given by Akin [Ref. 8]. 
The Type II boundary condition of equation (4.2) above is 
directly relatable to the fiber ends that move in unison, 
U, + (-1)U, = 0 


U, + (-1)U, =0 


& 


U, eae =140 (4.6) 


Une 1) U, = 0 


Zn 


tae a=-1 and c=0. The Type II boundary conditions of 
Equation (4.6) were applied in the stress analysis program. 
The final system of equations with both sets of boundary 
conditions applied was solved for the displacement vector U 
using the IMSL library subroutine DLSARG. 
Once nodal point displacements Us (j=1,2...,2n) were known 
fiber stresses go were determined using equation (2.14). This 


process was continued for each time step of output data 


Sys) 


obtained from the combustion code. The fiber stress analysis 


code is included in Appendix B. 
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V. CASE STUDIES 

The combustion code requires input data to start the 
analysis. Several input data files were produced encompassing 
various combinations of composite thickness, heat flux, 
differential pressure and combustion product fraction (F). 
Specific values of these variable parameters as used in the 
case studies are listed in Table 2. Several other important 
parameters were treated as constants for these studies as 
listed in Table 3. The input data files required by the 
combustion code were built covering all 1008 combinations of 
the variable parameters. These data files were processed by 
the combustion code and then the temperature profile data was 
further processed for stress failure information. Most of the 
lower heat flux results, primarily below 4000 btu/ eee 
produced relatively minor increases in plate temperatures 
which resulted in insignificant stresses and no failure. But 
the majority of runs gave interesting results. The final 
results of this effort are provided in graphical form in the 


following case studies. 
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TABLE 2: THE PARAMETERS VARIED IN THE CASE STUDIES. 


INPUT PARAMETERS THAT WERE VARIED FOR THE CASE STUDIES 
WITH THE VALUES INVESTIGATED 


COMPOS ITB. THICKNESS 2 Qi SM cei) cll ee 













CONSTANT HEAT FLUX 
(thousand btu/hr-ft?): 









DIFFERENTIAL PRESSURE (psi) 0270S 07> Ome 





COMBUSTION PRODUCT FRACTION (3%): 1) 0775.) 0 5 wom 





TABLE 3: PARAMETERS TREATED AS CONSTANT FOR THESE STUDIES. 





CONSTANT INPUT PARAMETERS TO THE COMBUSTION CODE 
FOR THE CASE STUDIES 


AMBIENT AIR TEMPERATURE: 80°F 
AMBIENT AIR PRESSURE: 2116.8 lb/ft? 
TORTUOSITY: 1.4 


PARTICLE SHAPE FACTOR: 272 


ORDER OF THE CHEMICAL REACTION: 0.5 


TIME STEP 
(between combustion code profiles): 0.5 seconds 


STOICHEOMETRIC RATIO: 02275 


eRe! 


A. TEMPERATURE PROFILES 

As stated in previously, induced stresses are dependent 
upon the fiber temperatures. Figures 11 through 16 show the 
temperature profiles of some representative combustion code 
runs up to the point of initial fiber failure. Note that the 
maximum temperature is always at the bottom of the plate where 
the heat flux is being applied and the air is coming in, while 
the cooler fibers are at the opposite surface. The main 
concept to bare in mind for these studies is that graphite 
fibers contract upon heating, so the hotter fibers will be in 
a tensile mode (see Figure 4). Since the tensile strength of 
these hot fibers is very high, these fibers will not fail. 
However, the fibers that are cooler will be in a compression 
mode with a low critical buckling strength and they will fail 
i buckling. 

A comparison of Figures 11, 12, and 13 show the effect of 
increasing the heat flux on a constant thickness material. 
Notice that for the lower heat fluxes, the temperatures are 
lower and the profile is more gradual. This, of course, is 
due to the lower heat flux applied but also due to the ability 
of the air flow to convectively transfer the heat away from 
the heated surface. If the heat flux is low, the flow of air 
through the medium has more time to distribute the heat 


through the material before failure occurs. If the heat flux 
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is high, the surface temperature rises quickly and heating is 
confined primarily to the area close to the heated surface. 

A comparison of Figures 11 and 16 indicates the changes 
in the temperature profile when the highest and lowest 
differential pressures are applied for a given heat flux and 
material thickness. As stated previously, the convective heat 
transfer in the medium has a large effect on the temperature 
profiles for the lower heat fluxes. For a constant thickness, 
the pressure gradient AP/L will be higher for a higher 
differential pressure. This higher pressure gradient causes 
a higher air flow through the material which increases the 
convective heat transfer. This increased heat transfer causes 
both a lower and more gradual temperature profile in the low 
heat flux cases. However, in the high heat flux cases this 
effect is insignificant as seen in Figures 12 and 14 for a 
heat flux ot 10,000nbEuy mete. 

As a final comparison Figures 12 and 15 are compared. 
These figures show the dependence of the temperature profile 
on thickness of the material while holding differential 
pressure constant. The temperatures at the surface of the 
1.0" plate are higher than the 0.5" plate. This can be 
explained as follows. If the AP is a fixed value, then the 
thicker plate has a smaller pressure gradient (AP/L) because 
the thickness L is larger. This smaller AP/L gradient will 
Cause less air flow through the thicker material resulting in 


a lower amount of convective heat transfer. The end result for 
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the thicker material is a higher and more steep temperature 
profile. This profile causes more fibers to be in the 


compressive mode and subsequently fail by buckling. 
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TEMPERATURE PROFILE 


WITH RESPECT TO TIME 
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Figure 11: Temperature profile up to first failure for 0.5", 
3,000 btu/hr ft* and 0.7 psi. 
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TEMPERATURE PROFILE 


WITH RESPECT TO TIME 


TIME TO FIRST FAILURE=12.1 sec 
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Figure 12: Temperature profile up to first failure for 0.5", 
morooo btu/hr ft* and 0.7 psi. 


43 


TEMPERATURE PROFILE 


WITH RESPECT TO TIME 
TIME TO INITAL FAILURE=4.9 sec 
HEAT FLUX=25,000 BTU/HR FT 
0.5 INCH THICKNESS 
AP=O.7PSI 
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Figure 13: Temperature profile up to first failure for 0.5", 
25,000 btu/hr ft? and 0.7 psi. 
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TEMPERATURE PROFILE 


WITH RESPECT TO TIME 


TIME TO FIRST FAILURE=12.0 sec 
HEAT FLUX=10,000 BTU/HR FT 
0.5 INCH THICKNESS 


AP=0.05PSI 
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Figure 14: Temperature profile up to first failure for 0.5", 
10,000 btu/hr ft* and 0.05 psi. 
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TEMPERATURE PROFILE 


WITH RESPECT TO TIME 


TIME TO FIRST FAILURE=27.3 sec 
HEAT FLUX=10,000 BTU/HR FT 
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Figure 15: Temperature profile up to first failure for 1.0", 
10,000 btu/hr ft* and 0.7 pou 
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TEMPERATURE PROFILE 


WITH RESPECT TO TIME 


TIME TO FIRST FAILURE=58.3 sec 
HEAT FLUX=3,000 BTU/HR FT 
0.5 INCH THICKNESS 
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Figure 16: Temperature profile up to first failure for 0.5", 
Beooo btu/hr ft* and 0.05 psi. 
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B. TIME TO INITIAL FAILURE (AP EFFECT) 

Figures 17 through 20 present the time to initial failure 
showing its dependency on both heat flux and differential 
pressure. As heat flux decreases, the time to initial failure 
progressively increases. For each specific AP there appears 
to be a lower bound of heat flux below which failure will 
never -occur- The figures also show that this lower bound 
value increases with increasing AP. However, this lower bound 
cannot be exactly determined through use of a computer code. 
Through the use of finite increments of heat flux, given 
previously in Table 2, this lower bound for failure was 
roughly determined. As seen in Figures 17 through 20 this 
lower bound of heat flux was found to within 500 btu/hr-ft?. 

For all thicknesses, the effect of pressure is relatively 
small except at the lower heat fluxes where the effect of air 
flow through the material has a large effect on the convective 
heat transfer. These pressure effects become ineffectual at 


higher heat fluxes. 
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TIME TO INITIAL FIBER FAILURE 
(EFFECT OF AP) 
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F=0.5 % 
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THIS HEAT FLUX 
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Figure 17: The time to initial failure (effect of 
differential pressure) for 1.0", and F=0.5%. 
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Figure 18: The time to initial failure (effect of 
differential pressure) for 0.75", and F=0.5%. 
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TIME TO INITIAL FIBER FAILURE 
(EFFECT OF AP) 
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Figure 19: The time to initial failure (effect of 
differential pressure) for 0.5", and F=0.5%. 
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Figure 20: The time to initial failure (effect of 
differential pressure) for 0.25", and F=0.5%. 
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C. TIME TO INITIAL FAILURE (F EFFECT) 

Figures 21 through 24 demonstrate how the critical 
buckling strength, and hence time to initial failure, varies 
with the combustion product factor (F). The value of F 
greatly effects not only the time to failure but determines 
whether or not failure will occur. As F increases there is 
more lateral support for the fibers so the critical buckling 
strength increases. This increase results in fewer fibers 
failing by buckling. In the case of the 0.25" plate none of 
the fibers reached the critical buckling strength when F=1.0%. 
All figures show that for a given heat flux, as F increases, 
the time to initial failure increases. 

The result of this analysis shows that the value assumed 
for the combustion product fraction (F) greatly influences the 
amount of time to initial failure of the material. The 
figures show that when F increases the resistance to buckling 


increases and buckling is less likely to occur. 
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Figure 21: The time to initial failure (effect of combustion 
product fraction) for 270" F"and AP=0.7 pee 
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Figure 22: The time to initial failure (effect of combustion 
product fraction) for 0.75", and AP=ji 7 spc 
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TIME TO INITIAL FIBER FAILURE 
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Figure 23: The time to initial failure (effect of combustion 
product fraction) for 0.5", and AP=0.7 psi. 
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Figure 24: The time to initial failure (effect of combustion 
Peoauce traction) for 0.25", and AP=0.7 psi. 
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D. PROGRESSION OF MATERIAL FAILURE (AP EFFECT) 

Figures 25 and 26 track the percentage of failed fibers 
from the time of the first observed fiber failures. Due to 
the 0.5 second time step used for temperature profile output, 
some of these curves start at time zero with a significant 
number failed fibers. Actually the failure represented by 
these curves commenced some time between the time of zero of 
the figures and the previous 0.5 second interval. Since 
observation of trends in this case are important, the 
increased computer time to determine a more precise starting 
time was not justified. 

The general trend is, the higher the differential 
pressure, the lower the percentage of failure. This is 
apparently is due to the more gradual temperature profiles 
caused by the higher air flows through the material. For the 
higher differential pressures/higher air flowrates, more 
fibers are heated to approximately the same temperature, and 
are therefore in the tensile mode prior to the time of initial 
fiber. These tensile fibers do not fail. 

Another item worth mentioning is the rapidity of failure. 
The materials progressed to their maximum percentage of fiber 
failure very rapidly. For the 1.0" thicker material this 
maximum failure percentage is reached within 3 seconds of the 
first fiber failure and for the 0.5" thinner material 
approximately 4 seconds. These times are independent of 


differential pressures. 
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Figure 25: The progression of material failure (the effect 
Oi Se) £or 0.5", and 5,000 btu/ hr-ft°. 
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Figure 26: The progression of material failure (the effect 
@EeAP)) for 1.0", and 5,000 btu/hr-ft?. 
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E. PROGRESSION OF MATERIAL FAILURE (EFFECT OF HEAT FLUX) 
Figures 27 and 28 indicate that as heat flux increases, 
time to maximum failure changes only marginally, effecting the 


3 to 4 second failure time by only 1 second. 
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PROGRESSION OF FAILURE 
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Figure 27: The progression of material failure (the effect 
Of: Neat flux) for 0.5", and AP=0.7 psi. 
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Figure 28: The progression of material failure (the effect 
Suemeat flux) for 1.0", and AP=0.7. 
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F. MAXIMUM FAILURE (EFFECT OF AP) 

Figures 29 through 32 are a consolidation the maximum 
fiber failure for material failure. These maximum percentages 
varied slightly with pressure at lower heat fluxes. As shown, 
for a given material thickness, higher heat fluxes give higher 
percentage of failed fibers. This is due to the more non- 
uniform temperature profile discussed previously. The maximum 
percentage failure seen in these studies is the case of the 
one inch plate at 25,000 btu/hr-ft?, which reached a maximum 


Of 97.56. 
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Figure 29: The maximum percentage of failure (the effect of 
me) for 1.0", and F=0.5%. 
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Figure 30: The maximum percentage of failure (the effect of 
mee) for 0.75", and F=0.5%.. 
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Figure 31: The maximum percentage of failure (the effect of 
AP) for. 0.5", and’) F=0. 5. 
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Figure 32: The maximum percentage of failure (the effect of 
AP) £Or 0.25", and. F-02524 
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G. MAXIMUM FAILURE (EFFECT OF F) 


The effect of the combustion product fraction is shown in 
Figures 33 and 34. As expected, assuming a higher value of 
the combustion product fraction (F) reduces the amount of 
failure because of the increase in buckling strength. For the 
0.25" plate there was no failure above a factor of 1.0% for 
range of heat fluxes investigated. 


These trends were independent of differential pressure. 
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Figure 33: The maximum percentage of failure (the effect of 
combustion product fraction) 0.5", and AP=0.7 psi. 
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Figure 34: The maximum percentage of failure (the effect of 
combustion product fraction) 0.25" andear=077 
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H. MAXIMUM FAILURE (EFFECT OF THICKNESS) 

As seen in Figures 35 and 36, the thickness of the 
material has minimal effect on the final percentage of 
material failure for higher heat fluxes. However at lower 
heat fluxes the thinner materials experience significantly 
less failure. In all cases the final percentage of fiber 


failure was between 42% and 97.5%. 
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Figure 35: The maximum percentage of failure (the effect of 
material thickness) F=0.5%, AP=0.7 psi. 
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Figure 36: The maximum percentage of failure (the effect of 
material thickness) F=0.5%, and AP=0.05. 
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VI. CONCLUSIONS 
The conclusions drawn from this study are as follows: 
The amount of time required before onset of fiber failure 
ranged from a very small 3 seconds (for a 0.25" plate at 
Pe go ey liebe) up to a maximum of 210 seconds (fora 
1.0" plate at 2,000 btu/hr-f£t*) . Given the reaction times 
of highly trained fire fighting teams onboard aircraft 
Carriers, it is quite conceivable that failure may be 
averted for the longer onset of failure times. 
Once the failure of carbon graphite-epoxy composite 
material begins maximum failure occurs very quickly in 
less than 4 seconds. This time is of such short duration 
that reducing the amount of maximum failure once failure 
begins is not possible. 
The predominant factors effecting survivability of carbon- 
graphite epoxy composite materials are the heat flux 
applied, the combustion product factor and the thickness 
of the material. The differential pressure effect varies 
from very important at low heat fluxes to insignificant at 
high heat fluxes. However, the effect of differential 
pressure is mitigated by the fact that response times are 
in all cases very small at the high heat fluxes. Once 
failure occurs, maximum failures in the range of 80-90% 


preclude any future use of the material. 
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APPENDIX A: COMBUSTION CODE EQUATIONS 


EQUATION SYMBOLS 


dynamic viscosity of air 
oxygen concentration 
air density 
graphite density 
specific heat capacity of air 
specific heat capacity of graphite 
Mass diffusion coefficient 
convection heat transfer coefficient 
air thermal conductivity 
effective thermal conductivity of porus graphite 
pseudo thermal conductivity to account for radiation 
pressure 
porosity 
oxygen reaction rate 
graphite reaction rate 
air temperature (°R) 
graphite temperature (°R) 
time 
pore velocity (of air) 


Specific internal area 
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FINAL FIELD EQUATIONS 


Following are the equations used in the Vatikiotis 
Combustion code. These equations are presented here in their 
form prior to non-dimensionalization. These equations were 
non-dimensionalized prior to use in the actual code. 


THE GRAPHITE HEAT TRANSFER EQUATION 


0 OT _ OT 
a, (oi-P) je aD i) ee - h2(T, Sal Ale 1K a = (1-p) Pg ae 


THE INTERNAL AIR HEAT TRANSFER EQUATION 

















0 oT, or, od oft Sy O(pP) _ OT, 
x (Pka Gy) ~ PPagCat ay * BACT" Te) * UWS = Pala ae 
OXYGEN TRANSPORT EQUATION 
96 » d > _ _, 
= (PP ae) Ox en eH o2i7 PSt 
THE INTERNAL FLOW EQUATION 
gee. pL SPe By | 2 GM) ela) Gee, 
dx? p, Ox m Ox Hb Ox dx Saino G 
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BOUNDARY CONDITIONS IMPOSED 


For the X=0 location, 


T, = T. 
> = 6. 
P= P, 

During neat aqqi Cio \ b-o) cae) sie = Sen 


After fire is removed (k, + K,) = aaa Pee oi) eas oe(T, — fap 


For the X=L location, 





OT 
(Ko KD Se = yt, Te) Coes 
OT, =o 
Ox 
Op _ 
aeiaee 


where: 
Go, = the heat £lusc apeilied 
T. = ambient temperature (°R) 
d?,. = ambient oxygen concentration 


h, and hs = convection heat transfer coefficients 


68 


INITIAL CONDITIONS 
There are a number of ways to start these "combustion" 
problems. One way is to set the oxygen concentrations and the 
air and graphite temperatures to ambient conditions then apply 
an external heat source. This "natural" starting scheme 
requires some amount of computer effort to bring the plate up 
to conditions where temperatures have increased significantly. 
Instead of starting the problem at ambient conditions, 
computational effort can be reduced by starting the problem 
with the plate at elevated temperatures. 
In this particular study the initial conditions were; 
Graphite temperature 400°F 
Air temperature 400°F 


Oxygen concentration 0.0172 lb/ft? (ambient level). 
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APPENDIX B: STRESS ANALYSIS CODE 


*EDWARD FAXLANGER JULY 1992 


wk we Ke Ke WK KKK HK KK KKK KKK KI KKK KKH HK HHH KKK KK KKK KEKEKKKK KKK KKK KK KKK KKK KK 


* 
* 


THIS PROGRAM CALCULATES STRESSES INDUCED IN LONGITUDINAL FIBERS 
RESULTING FROM TEMPERATURE VARIATION ACROSS THE MEDIUM 


KKK KKkKkKKKKKKKKK KKK KKK K KKK KK KKK KKKKKKKK KKK KK KKK KKK KKK KKK KK 
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DESCRIPTION OF VARIABLES 


U = DISPLACEMENTS OF NODAL POINTS 
BC = VALUE OF THE TYPE I BOUNDARY CONDITION FOR NODAL POINT IBC 
A = MATRIX CONTAINING INDIVIDUAL STIFFNESS 
B = MATRIX CONTAINING INDIVIDUAL FORCE VECTOR 
AA = MATRIX CONTAINING GLOBAL STIFFNESS 
BB = MATRIX CONTAINING GLOBAL FORCE VECTOR 
Y = LOCATIONS OF THE INDIVIDUAL FIBER (FROM BOTTOM OF PLATE) 
TEMP = INTERPOLATED TEMPERATURE OF THE FIBER 
DIAM = INTERPOLATED DIAMETER OF THE FIBER 
SIGMA = CALCULATED STRESS OF THE INDIVIDUAL FIBER 
YIN = LOCATIONS OF NODAL POINTS FROM THE COMBUSTION CODE 


TEMPIN=TEMPERATURES FROM THE COMBUSTION CODE 
DIAMIN=THE DIAMETER AT THE NODAL POINT FROM THE COMBUSTION CODE 


IBC = THE NUMBER OF THE NODE HAVING A TYPE I BOUNDARY CONDITION 
NODES= A MATRIX USED FOR ELEMENT CONNECTIVITY TO THEIR NODAL POINTS 
Yl = LOCATION OF THE NODAL POINTS BROUGHT IN FROM VATIKIOTIS CODE 
NFLAG = A FLAG TO INDICATE TYPE OF FAILURE 


NELEM =A MATRIX CONTAINING THE NAMES OF THE UNBROKEN FIBERS 
NE = NUMBER OF FIBER ELEMENTS 

NNODE = TOTAL NUMBER OF FIBER NODES 

PERCNT= PERCENTAGE OF ORIGINAL BUCKLING STRENGTH USED 

TSTREN= TENSILE STRENGTH OF FIBERS 


CSTREN= COMPRESSIVE STRENGTH OF FIBERS (NO BUCKLING) 
EFIB = ELASTIC MODULUS OF THE FIBERS 
VF = FIBER VOLUME FRACTION (ASSUMED CONSTANT) 
EMAT = ELASTIC MODULUS OF THE MATRIX 
POISM = POISON RATIO OF THE MATRIX 
ALPHA = THERMAL EXPANSION COEFFICIENT (AXIAL) 
WIDTH = THE ASSUMED LENGTH OF EACH FIBER (ASSUMED TO BE ONE) 
TEMREF= REFERENCE TEMPERATURE (WHERE STRAIN EQUAL ZERO) 
NNP = THE NUMBER OF NODAL POINTS READ IN FROM COMBUSTION CODE 
TIME = TIME OF THE TEMPERATURE PROFILE BEING ANALYZED 
THICK = THE ORIGINAL THICKNESS OF THE MATERIAL (BEFORE BURNING) 
HEIGHT= THE THICKNESS OF THE MATERIAL AT THE PRESENT TIME 
NEUSED= THE NUMBER OF UNBROKEN FIBERS AT PRESENT TIME 
NELEM = A STORAGE OF UNBROKEN FIBERS NAMES 
GMAT = SHEAR MODULUS OF THE EPOXY MATRIX MATERIAL 


MAIN PROGRAM 
IMPLICIT REAL*8 (A-H,0-Z) 
THE FOLLOWING PARAMETER N1 SETS THE DIMENSION OF ARRAYS AND 
NEEDS TO BE EQUAL TO THE HIGHER OF 1)THE NUMBER OF NODES BEING 
READ IN **OR** 2) THE NUMBER OF ELEMENTS BEING ANALYZED 
PARAMETER (N1=101) 
PARAMETER (NN=2*N1) 
COMMON STATEMENT FOR MATERIAL PROPERTIES 
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COMMON / PROPER/TSTREN, CSTREN, EFIB, EMAT, ALPHA, WIDTH, THICK, POISM, VF 
& 
REAL*8 U(NN) ,BC(NN) ,A(2,2),B(2),BB(NN) ,AA(NN,NN) ,Y(N1), 
&TEMP (N1) ,DIAM(N1) ,SIGMA(N1) , YIN(N1) , TEMPIN(N1) , DIAMIN(N1) 
DIMENSION IBC (NN) , NODES (NN, NN) , INDEX (2) , NFLAG(N1) , NELEM (N1) 


*THE FOLLOWING 3 STATEMENTS INCREASE THE NUMBER OF EQUATIONS 
*THAT THE EQUATION SOLVER DLSARG WILL HANDLE. IF LARGER SPACE IS 
* NEEDED DLSARG WILL GIVE AN ERROR MESSAGE AND RECOMMEND LARGER 
* WORK SPACES 

COMMON /WORKSP/RWKSP 

REAL RWKSP (82236) 

CALL IWKIN (82236) 
C DATA FILE FOR OUTPUT OF DATA 

OPEN (UNIT=11, FILE='’/TEST31 STRESS B1’ , STATUS=' OLD’ ) 
C DATA FILE FOR READING INPUT DATA 

OPEN (UNIT=12, FILE=’ /TEST31 TEMPER Bl’ ,STATUS=’ OLD’ ) 
C NE IS THE NUMBER OF FIBERS ASSUMED FOR ANALYSIS 

NE=30 

NNODE=NE* 2 
C ENTER THE PERCENTAGE OF MATRIX SUPPORT 

PERCNT=.005 


C FIBER CHARACTERISTICS; LONGITUDINAL COMPRESSIVE AND TENSILE 
C STRENGTHS (PSI), ELASTIC MODULUS AND THERMAL COEFFICIENT 
TSTREN=300 .0D3 
CSTREN=-260.0D3 
EPs =sd),0D6 
EMAT=.5D6 
VF=. 196 
POLSM= 35 
ALPHA=-0.55D-6 
C WIDTH OF SIMULATED HOLE (LENGTH OF FIBERS) ASSUMED TO BE UNITY 
C SINCE IT HAS NO EFFECT ON PROBLEM 
WIDTH=1 .0DO 
C REFERENCE TEMPERATURE IN DEGREES F WHERE INITIAL STRESSES WERE ZERO 
TEMREF= 60.D0 


C INITIALIZE FLAGS TO ZERO (NFLAG KEEPS TRACK OF BROKEN FIBERS) 
e (NSTOP STOPS THE PROGRAM ) 
NSTOP=0 
DO M=1,NE 
NF LAG (M) =0 
END DO 


* READ IN OF INFORMATION FROM COMBUSTION CODE OUTPUT 
READ (12, *) NNP 
CALL CHECK (NNP, NSTOP) 
IF (NSTOP.NE.0)GOTO 100 
CALL READER (N1,NNP, TIME, HEIGHT, YIN, TEMPIN, DIAMIN, NSTOP) 
IF (NSTOP.NE.0)GOTO 100 


*LOOP FOR FILLING ALL ELEMENT HEIGHTS ASSUMING EVEN SPACING 
THICK=HEIGHT 
YINC=THICK/ (NE) 
YSTART=YINC/2 
DO M=1,NE 
Y (M) =YSTART+YINC* (M-1) 
END DO 


qo) 


200 BURN=THICK-HEIGHT 
C INITIALIZE ARRAYS TO ZERO 


DO f= hnNN 

BB (I) =0.0D0 
DOI= NN 

AA (IT, J) =0, 0D0 
END DO 

END DO 

DO I=1,NE 


C IF THE FIBERS ARE ALREADY CONSUMED BY COMBUSTION (OUTSIDE EDGE) 
C NULL THEIR EFFECT BY SETTING THEIR FLAG TO 1 
IF (Y(1I) .LT.BURN) THEN 
NFLAG (I) =1 
END IF 
DETERMINE THE TEMPERATURE AND DIAMETER OF NON-FAILED FIBERS THIS 
LOOP IS FOR INTERPOLATION OF DATA RECEIVED FROM "STRESS INPUT" FILE. 
THE RECEIVED DATA IS IN RELATIVE SPACING AND MUST BE CHANGED TO 
ABSOLUTE SPACING IN ORDER TO DETERMINE PERCENTAGE FAILURE OF FIBERS 
INTERPOLATE NON-BROKEN FIBERS ONLY 
IF (NFLAG(I) .EQ.0) THEN 
DO N=1,NNP-2 
IF((Y(I) .GE.YIN(N) +BURN) .AND. (Y(I) .LE.YIN(N+2) +BURN) ) THEN 
CALL INTERP(YIN(N) +BURN, YIN (N+1) +BURN, YIN(N+2) +BURN, 
&DIAMIN(N) , DIAMIN (N+1) , DIAMIN(N+2) , Y(I) , DIAM(I) ) 
CALL INTERP(YIN(N) +BURN, YIN (N+1) +BURN, YIN (N+2) +BURN, 
&TEMPIN(N) , TEMPIN(N+1) , TEMPIN(N+2) , ¥ (I) , TEMP (I) ) 
END IF 
END DO 
END IF 
END DO 


Or GO 


C LOOP FOR DETERMINING THE ELEMENTS THAT ARE NOT BROKEN (IE THE 
C ONES THAT ARE TO BE ANALIZED) 

NEUSED=0 

DO I=1,NE 

IF (NFLAG (I) .EQ.0) THEN 

NEUSED=NEUSED+1 

NE LEM (NEUSED) =I 

END IF 

END DO 

IF (NEUSED .LT.2) THEN 

WRITE (11,*) ’ALL BUT TWO ELEMENTS HAVE BROKEN’ 

WRITE (11,*) ‘ANALYSIS IS STOPPED/PROGRAM STOPPED’ 

GOTO 100 

END IF 

NNUSED=NEUSED* 2 


* LOOP FOR ESTABLISHING CONNECTIVITY OF NODES TO THEIR ELEMENTS 
* FOR USE IN BUILDING GLOBAL MATRIX AND GLOBAL VECTOR 

L=1 

DO I=1,NEUSED 

DO J=1,2 

NODES (I,J) =L 

L=L+1 

END DO 

END DO 


* LOOP TO BUILD GLOBAL MATRIX (AA) AND VECTOR (BB) 
DO IE=1,NEUSED 
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CALL MATRIX (A,B, TEMP (NELEM(IE) ) -TEMREF, DIAM (NELEM (IE) ) ) 
INDEX (1) =NODES (IE, 1) 

INDEX (2) =NODES (IE, 2) 

DO I=1,2 

ITI =INDEX (I) 

BB (II) =BB (II) +B (TI) 

BIO) iil % 

JJ=INDEX (J) 
AA(II,JdJ) =AA(II, JJ) +A(I,J) 
END DO 

END DO 

END DO 


*LOOP FOR TYPE I BOUNDARY CONDITION INITIALIZATION 
PO r=1 NNUSED-i, 2 
1) aE We) 
BC (I) =0.0D0 
END DO 


* APPLY TYPE 1 BOUNDARY CONDITIONS 
DO I=1,NNUSED-1, 2 
USIE EIS (0), 

DO K=1,NNUSED 

BEM) BE (K) = BC (1)i*AA(K, J) 
AA (K,d) =0.0 

AA(J,K) =0.0 

AA(J,d) =1.0 

BEng) Sein) 

END DO 

END DO 


C APPLY TYPE 2 BOUNDARY CONDITIONS 
Desd—2 NNUSED-2; 2 
CALL MOD2FL(NN,NNUSED,J,J+2,-1.0D0,0.0D0,AA, BB) 
END DO 


C SOLVE THE SYSTEM OF EQUATIONS USING IMSL LIBRARY 
CALL DLSARG (NNUSED, AA, NN,BB,1,0U) 


DO I=1,NE 

SIGMA (I) =0.0D0O 

END DO 

DO I=1,NEUSED 

STRAIN=U (I*2) /WIDTH 

SIGMA (NELEM (I) ) =EFIB* (STRAIN-ALPHA* (TEMP (NELEM(I) ) -TEMREF) ) 
END DO 


C LOOP TO CALCULATE STRESSES OF EACH ELEMENT AND OUTPUT RESULTS 
ee20Q OUTPUT FILE 

WRIDE (iio) Tove 

WROTE (11776) 

DO I=1,NE 

IF (NFLAG (I) .NE.0)WRITE(11,8) I,NFLAG(I) 

IF (NFLAG (I) .EQ.0)WRITE (11,7) I,SIGMA(I) , DIAM(I) 

END DO 

FORMAT (SX, ‘FIBER NUMBER’ ,5X,’'FIBER STRESS (PSI)‘’,2X, ‘DIAMETER’ ) 
PORMAT (2 logit, GlOe4s)6x,G10.4,3X,G11.4,2X,11) 

FORMAT (17X,13,4X, ‘PREVIOUSLY BROKEN, NFLAG=’ ,1I1) 

BORMATR OL vwEMb= F823, SECONDS’ ) 

CALL BROKEN (N1,NE, SIGMA, DIAM, NFLAG, BROKE, PERCNT) 


WRITE (11, *) 1% 2% ee ee Fe hee te ie ek II I II IOI II II IO I Ie dee te de def 


(9 © ~) OF 
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WRITE (11733 
CALL READER (N1,NNP, TIME, HEIGHT, YIN, TEMPIN, DIAMIN, NSTOP) 
IF (NSTOP.NE.0)GOTO 100 
GOTO 200 
100 END 


Cc 
Cc 
SUBROUTINE MATRIX (A,B, TEMP, DIAM) 
C THIS SUBROUTINE CALCULATES THE ELEMENT LEVEL 
C STIFFNESS MATRIX (A) AND VECTOR (B) 
IMPLICIT REAL*8 (A-H,0-Z) 
COMMON /PROPER/TSTREN,CSTREN,EFIB, EMAT, ALPHA, WIDTH, THICK, POISM, VF 
REAL*8 A(2,2),B(2) 
PI=4*ATAN (1.0) 
AREA=DIAM*DIAM*PI/4.0 
S=AREA*EFIB/WIDTH 
A(1,1)<=S 
A(1,2)=-S 
A(2,1)=-S 
A(2,2)2S 
B (2) =EFIB*AREA*ALPHA* TEMP 
BCL) - Biz) 
END 


Qa 


SUBROUTINE BROKEN (N1,NE,STRESS, DIAM, NFLAG, BROKE, PERCNT) 
THIS SUBROUTINE DETERMINES THE NUMBER OF BROKEN FIBERS (BROKE) 
DUE TO TENSILE AND COMPOSITE BUCKLING 
NFLAG EQUALS 0:FOR UNFAILED FIBER 

1:FOR BURNT AWAY FIBERS 
2:FOR OVER STRESS BY TENSION 
3:FOR OVER STRESS BY BUCKLING 

IMPLICIT REAL*8 (A-H,0-Z) 

COMMON /PROPER/TSTREN, CSTREN, EFIB, EMAT, ALPHA, WIDTH, THICK, POISM, VF 

REAL*8 STRESS (N1) , DIAM(N1) 

DIMENSION NFLAG (N1) 

PI=4*ATAN (1.0) 

C CALCULATE COMPOSITE BUCKLING STRENGTH OF THE FIBER 
C AS A FRACTION OF THE ORIGINAL COMPOSITE BUCKLING STRENGTH 
GMAT=EMAT/ (2.* (1.+POISM) ) 
BUCKLE=-GMAT/ (1.-VF) *PERCNT 
DO I=1,NE 
C DO NOT CHECK PREVIOUSLY BROKEN FIBERS 
IF (NFLAG(I) .EQ.0) THEN 
C CHECK FOR TENSILE STRESS FAILURE 
IF (STRESS (I) .GE.TSTREN) THEN 
NFLAG (I) =2 
END IF 
C CHECK FOR BUCKLING FAILURE 
IF (STRESS (I) .LE.BUCKLE) THEN 
NFLAG (I) =3 
END IF 
END IF 
END DO 
C SUMMATION OF EACH TYPE OF BREAKAGE 
BURNT=0.0 
TENSIL=0.0 
BUCK=0 .0 
DO I=1,NE 
IF (NFLAG (I) .EQ.1) BURNT=BURNT+1 
IF (NFLAG(I) .EQ.2) TENSIL=TENSIL+1 


oe Eis i @ lal @ 1 @ 
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CO er) 


ele ae 


212 Le ie he Ae ie ae a ee eo) 


QAM 


IF (NFLAG(I) .EQ.3) BUCK=BUCK+1 
END DO 
TBROKE=BURNT+TENSIL+BUCK 
WRITE (11,5) BURNT*100/NE 

5 FORMAT (1X, 'PERCENTAGE BURNT FIBERS’ ,F5.1) 
WRITE (11,4) TENSIL*100/NE 

4 FORMAT (1X, ‘PERCENTAGE BROKEN IN TENSILE MODE’ ,F5.1) 
WRITE (11,2) BUCK*100/NE 

2 FORMAT (1X, '’PERCENTAGE OF BUCKLED FIBERS’ ,F5.1) 
WRITE (11,1) TBROKE*100/NE 

1 FORMAT (1X, 'TOTAL PERCENTAGE OF FAILED FIBERS’ ,FS5.1) 
END 


SUBROUTINE INTERP (X1,X2,X3,Y1,Y2,Y3,XX, YY) 
THIS SUBROUTINE TAKES 3 X VALUES, 3 Y VALUES. THEN USING A GIVE XxX 
VALUE PERFORMS QUADRATIC INTERPOLATION TO DETERMINE THE YY VALUE. 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 A(3,3),B(3),C(3) 

PA 1) =12 D0 

A(2,1)=1.D0 

Als) 2120 


A(2,2) =X2 
A(3,2) =X3 


A(1,3) =X1**2. 
7 PARE) EP PA oa 
A(3,3) =X3**2. 


B(1)=Y1 
B (2) 
B(3) = 73 


CALLING IMSL LIBRARY SUBROUTINE TO SOLVE A(3,3) *C(3) =B(3) 
TO DETERMINE THE QUADRATIC COEFFICIENTS 

CALL DLSARG(3,A,3,B,1,C) 

VY =€() G2) * xX + (3) xk * 27 

END 


SUBROUTINE MOD2FL(NN,NSIZE,L1,L2,A,B,S,C) 
IMPLICIT REAL*8 (A-H,0O-Z) 
DIMENSION S(NN,NN), C (NN) 
ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee 
APPLY TYPE 2 CONSTRAINTS TO FULL EQUATIONS 
METHOD AND SUBROUTINE BY: J. ED AKIN 
ee ee ee ee ee ee ee 
NN = THE TOTAL SIZE OF THE PASSED 
NSIZE = WORKING SIZE OF ARRAYS 
S = FULL SQUARE MATRIX 
C = COLUMN MATRIX 
L1,L2 = CONSTRAINED DEGREES OF FREEDOM SUBCRIPTS OF D MATIRX 
A,B = CONSTRAINT COEFFICIENTS 
Diets A ile). = Bo. S*D =C 
INITIAL CALCULATIONS 
Siga Ss (lah 12 ) 
Clia= CULT) 
SUBTRACT A*COLUMN L1 FROM COLUMN L2 
SUBTRACT B*COLUMN L1 FROM C 
DO 10 I = 1,NSIZE 


S (Glia?) Ss S (L712) 7 A*S (I, L1) 
Gi) =2e (t= - Bes (ti) 


75 


Gl @ MT @ i 


10 


Se) Oro 

Sueph IG) Se 

CONTINUE 

SLi le 

§$(L2,L2) = §(L2,L2) - A*S12 + A*A 
S42 eee 

S( ily bo) ee 

C{Lil} = 5 

C(L2) = C(bL2)" = A*crites. 
RETURN 

END 


SUBROUTINE READER (N1,NNP,TIME,HEIGHT, YIN, TEMPIN, DIAMIN, NSTOP) 


THIS SUBROUTINE READS IN THE NODAL POINT VALUES FROM THE 
COMBUSTION CODE 


IMPLICIT REAL*8 (A-H,0-2Z) 

REAL*8 YIN(N1) , TEMPIN(N1) , DIAMIN(N1) 
READ (12, *, END=10) 

READ (12, *, END=10) 

READ (12, *, END=10) 

READ (12,*,END=10) TIME 

READ (12, *, END=10) 

READ (12,*,END=10) HEIGHT 

READ (12, *, END=10) 


CONVERT HEIGHT TO INCHES 


HEIGHT =zBETGHT*12 20 
DO I=1,NNP 
READ (12,*,END=10) K,YIN(I),TEMPIN(I) , DIAMIN (I) 


CONVERT INPUT LOCATION(YIN) AND DIAMETER(DIAMIN) TO INCHES 


10 
20 


YIN(I) =YIN(I) *12.0 
DIAMIN (I) =DIAMIN (I) *12.0 
END DO 

GOTO 20 

NSTOP=1 

END 


SUBROUTINE CHECK (NNP, NSTOP) 


QUICK CHECK TO SEE IF THE NUMBER OF INPUT NODAL POINTS IS ODD 
SO QUADRATIC INTERPOLATION CAN BE USED. IF NOT STOP PROGRAM 


MCHECK=MOD (NNP, 2) 

IF (MCHECK.EQ.0) THEN 

WRITE (*,2)’THE NUMBER OF INPUT NODAL POINTS MUST BE AN ODD’ 
WRITE (*,2) NUMBER IN ORDER FOR THIS PROGRAM TO USE QUADRATIC’ 
WRITE (*,2) ‘INTERPOLATION. CHECK NNP IN INPUT DATA’ 

WRITE (*,2)’THIS PROGRAM IS STOPPED’ 

FORMAT (1X, A) 

NSTOP=1 

END IF 

END 
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APPENDIX C: SAMPLE OUTPUT DATA FROM STRESS CODE 


AbD IDS 0.000 SECONDS 


FIBER NUMBER FIBER 


WODTINUM FP WNPH 


ra 
O7 
OooueqocT0eoqcoqooqcooqoq0oqCcoeoooqooqoo 6 ooo 


PERCENTAGE BURNT FIBERS 0 


SVS SS) 


0 


PERCENTAGE BROKEN IN TENSILE MODE 


PERCENTAGE OF BUCKLED FIBERS 


TOTAL PERCENTAGE OF FAILED 


kakkkkkkkekkkkkkkkkkkkkkekekkekekkkekekekkekekekekkkkk 


0.0 
FIBERS 


(PST) 
.0000E+00 
.0000E+00 
.0000E+00 
.0O000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
. OO00E+00 
.0000E+00 
.0000E+00 
.0OO000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 
.0000E+00 


0.0 


0.0 


Ua 


DIAMETER 


Oo0oo0oocqooqoqoqocoqo fo 0 @C @ @@ Oe @ ee © eo eG ea ea & 


218) 1 CYS OP 
Hoes =02 
SooSom- 02 
moos E-O02 
Boos sh - O02 
Peso 38-02 
nooo 58-02 
peo som 02 
BeOS om U2 
moos sn- 02 
neo son O02 
ee SS 56> U2 
Hososn-> 02 
2353558 -02 
~O5338-02 
woo 508- 02 
.8333E-02 
O33 58-02 
.8333E-02 
.8333E-02 
Poss 5n—-OZ2 
Boe sor 02 
Pos53s5-02 
ooo oH -02 
B1e)6) Si Gls = 10% 
~6333E-02 
78S 558- 02 
HoSs5n-O2 
FOS S55 U2 
poo S 3h - OZ 


TIME= &. 310 = SheGOnps 


FIBER NUMBER FIBER STRESS 
Jl Oe Dae Oo 
2 2es0. 
a Loi. 
4 dbo ys hoy 
5 OZ Ze 
6 70S 
7 433.8 
8 206.0 
g Sao 7 

LO ss 1.2 
sae =249..3 
ZZ -340.1 
i -408.6 
14 =4259. 2 
i> ao). 9 
EG =e ee 2 eC) 
Ig 540.2 
18 2552.6 
19 =561..0 
20 =566.6 
Zed somo. 3 
Page -572.6 
26 e574". 1 
24 StL 
20 i. 7 
26 =576;..0 
Za =o) OS 
28 -576.4 
29 = S16. 5 
30) = 5 OF..9 


PERCENTAGE BURNT FIBERS 0.0 

PERCENTAGE BROKEN IN TENSILE MODE 
PERCENTAGE OF BUCKLED FIBERS 0.0 
TOTAL PERCENTAGE OF FAILED FIBERS 


wkKkkkkkkkkk kkk kkk kkk kK KKK KKK KKK KKK KKK KK KK 


(PST) 


O20 


O20 


Tee! 


DIAMETER 


oooooococooeo°ocooeoodoe0odococo0ocooeoeeceoc0°000o0ooo000o00 0c fo 


~8333E-02 
Tosser -O2 
.8333E-02 
763953 E02 
,8333E-02 
.8333E-02 
8333-02 
Foss om—O2 
~os54h- 02 
-o333h-02 
POSS Sma 
. 8333E-02 
.8333E-02 
733 Sez 
Poso3SE-O2 
24333E-02 
Poss 02 
Pos SSh-=02 
70545 5E 702 
FOSS SE 02 
SoosS5h—-o2 
.8333E-02 
Foss SH OZ 
Foss SE 02 
Foss Se - O02 
13334 -02 
FooSen- UZ 
783338-02 
pOS53H- 02 
Toso S ne OZ 


MiME =e? 610 SECONDS 


FIBER NUMBER FIBER STRESS 


1 3524. 
2 BO22. 
3 ae) 
4 Zo 0:: 
iS) £3805. 
6 1443. 
7 Oe: 
8 797.38 
2 515.9 
10 ZGay 7 
ital 34.97 
Wee = Anya ae 
13 =O 
14 -490.6 
ZS -619 24 
16 eee Orr 
alg) foals Bea! 
18 -894 .3 
9 2955.8 
20 =1006- 
21 -1046. 
22 ee 7s 
23 = IL SO ee 
24 LZ. 
25 eel 3 Gre 
26 -1148. 
Pas = HALES « 
28 = akalisyal 
20 ald Go 
30 wee 7 


PERCENTAGE BURNT FIBERS 0.0 
PERCENTAGE BROKEN IN TENSILE MODE 
PERCENTAGE OF BUCKLED FIBERS 13.3 


(PST) 


0 


0) 


TOTAL PERCENTAGE OF FAILED FIBERS 13.3 


wkkekkkkkkkkkkkkkkkk kk kkk kk kkk kk kk kk kk kkk kk 


19 


DIAMETER 


CPoo0o00e0o0eoecoce0cococoeco0c0co0o0ooo°0oeq0oo0oo0oc ooo coca 080 


Pao ssh OZ 
Seas Sr- O02 
moo ooE-OZ 
oS Son-O2 
moss 5h -O02 
BOS oor-O02 
Boo S3E-02 
FOS S35E-O02 
HOS 53E-02 
moo SSE 02 
mos 55-02 
eoo 356-02 
oS SOs OZ 
mes338-02 
Poo somo 2 
meso sh -O2 
oes 5h-O2 
7O5S3E-02 
283538 -02 
Zoosk -O2 
(eS 2) 
woos Se-02 
eos558-02 
moss sb-02 
weoSsk-O02 
mes33h-O02 
Boo o on -O2 
.8333E-02 
26333E-02 
Seo ss OZ 


TIME= 26.560 SECONDS 

FIBER NUMBER FIBER STRESS 
i aS 
Z eo: 
3 1210). 
4 ieo2.. 
5 1174. 
6 902:0 
i SO0SieZ 
8 29755 
9 =S7231 
10 2290.6 
a Gal =561.4 
12 =O be o 
i -1044. 
14 SoS. 
LS -1445. 
16 TALS shiere 
7 -1768. 

18 PREVIOUSLY 

i PREVIOUSLY 

20 PREVIOUSLY 

Zu PREVIOUSLY 

PAP) PREVIOUSLY 

as PREVIOUSLY 

24 PREVIOUSLY 

25 PREVIOUSLY 

26 PREVIOUSLY 

Pad PREVIOUSLY 

28 PREVIOUSLY 

29 PREVIOUSLY 

30 PREVIOUSLY 

PERCENTAGE BURNT FIBERS 0.0 


PERCENTAGE BROKEN IN TENSILE MODE 
PERCENTAGE OF BUCKLED FIBERS 56.7 
TOTAL PERCENTAGE OF FAILED FIBERS 


Ke te eke te te te te te te te te te te te ke oe te oie te te ie keke ete te keke te ete ke keke ke eke ee te 


(PSY) 


BROKEN, 
BROKEN, 
BROKEN, 
BROKEN, 
BROKEN, 
BROKEN, 
BROKEN, 
BROKEN, 
BROKEN, 
BROKEN, 
BROKEN, 
BROKEN, 
BROKEN, 


Oo 


S67 
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DIAMETER 


Oooqoqceoo 0 OO Oo eo ee Ses 


20335 B02 
7O3SSSE02 
78333 -02 
-o34338-02 
Pos S3Eao72 
283328 -02 
~oSoSol- 02 
28333E-02 
-ooSs3h U2 
oo SSm 02 
FOS S33 C2 
POs 3on-o2 
gto CLEISYS GG 
FO353E 02 
7o5 3536-02 
TOS SSE -U¢ 
nosso E v2 


NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 
NFLAG=3 


APPENDIX D: SAMPLE INPUT DATA TO COMBUSTION CODE 


80.E0 

50.4E0 

21 16.8E0 
1.4E0 


Wees66055555555555507E-03 
so ciI0 5550555555550 7E-03 
0.416666666666666667E-01 
0.416666666666666667E-01 


53 .34E0 
86 .0E0 

Oe 231E0 
70. 3E0 

O02 9E0 
77 LEO 
Hr OEO 
14085 
0.5E0 
0.5E0 
07.37 5E0 
2.065E6 


28836 .E0 
30000.E0 


4.D0 
0.E0 
200 .E0 
40 .E0 

i 00OD-5 
ee 10 
Oe 50/0 
0.0010 
2500.E0 


.EO 


TAMB AMBIENT TEMPERATURE 
UINF DIFF. PRESS. ACROSS MEDIUM LB/FT2 OR SURF VELOCITY 
PAMB AMBIENT PRESSURE LB/FT2 
TORT TORTUOGS@ETY OF THE MEDIUM 
DINT LUMP FIBER DIAMETER PER CELL 
SINIT DISTANCE BETWEEN LUMPED FIBERS 


XLENGTH TOTAL MEDIUM THICKNESS 
XINIT INITIAL THICKNESS OF THE MEDIUM 


RGAS THE GAS CONSTANT FOR AIR 

GCOND THERMAL CONDUCTIVITY OF GRAPHITE 

SPHTG SPECIFIC HEAT OF GRAPHITE 

ROG DENSITY OF GRAPHITE 

EMIS THE EMISSIVITY OF THE GRAPHITE 

PHI SHAPE FACTOR 1.0 FOR SPHERE .91 FOR CYLINDERS 
HLENTH REF LENGTH FOR SURFACE HEAT TRANSFER (NORMALLY 1) 
ENTHAL THE ENTHALPY OF CHEMICAL REACTION 

ORDNTH THE ORDER OF THE CHEMICAL REACTION 

EXTEMP THE EXPONENT OF THE TEMPERATURE FOR REACTION RATE 
STIOCH STOICHIOMETRIC RATIO FOR THE REACTION 

RECOEFF ARRHENIUS EQUATION COEFFICIENT 

ENERGY ACTIVATION ENERGY DIVIDED BY THE IDEAL GAS CONSTANT 
SQ THE HEAT FLUX (SIMULATING THE FIRE) 

7 THE APPROXIMATE TIME INTERVAL FOR PRINTING RESULTS 
TSTART THE STARTING TIME OF THE PROBLEM 

TSTOP THE STOP TIME OF THE PROBLEM 

TQ THE HEATER TURN OFF TIME (FIRE IS OUT) 

HSTART THE INITIAL TIME INTERVAL IN SUBROUTINE SDESOL 
HMIN THE MINIMUM TIME INTERVAL FOR SUBROUTINE SDESOL 
HMAX THE MAXIMUM TIME INTERVAL FOR SUBROUTINE SDESOL 
RMSEPS ERROR CRITERIA FOR SUBROUTINE SDESOL 

TMAX THE MAX TEMP ON PRINTER PLOTS (Y AXIS) 

TMIN THE MIN TEMP ON PRINTER PLOTS (Y AXIS) 

CMAX THE MAX OXYGEN CONC. ON PRINTER PLOTS (Y AXIS) 
CMIN THE MIN OXYGEN CONC. ON PRINTER PLOTS (YAXIS) 

SL1 SLOPE FOR A LINEAR PROFILE FOR GRAPHITE TEMP 

SL2 SLOPE FOR A LINEAR PROFILE FOR AIR TEMP 

SL3 SLOPE FOR A LINEAR PROFILE FOR CONCENTRATION 


MAX ORDER EQUATION SOLVER USED IN SDESOL 


PRINT 


PARAMETER FOR SDESOL 


800 .E0 


OCP FODODOONHFONMNOCOO: . 
ye) e J @ e 
4 Hi 
29) es) 
Q C 
63 
Kc 


IQ 


MULTIPLIER OF PRINT TIME STEP ONCE PROGRAM IS IN SURF RECESS 
NUMBER OF NODAL POINTS 

O=AUTO 1=USER INPUT(SEE BELOW) OF NODAL LOCATIONS 

O=AUTO 1=USER INPUT(SEE BELOW) OF PARTICLE DIAMETER 


O=UINF IS DIFFERENTIAL PRESSURE 1=UNIF IS SURFACE VELOCITY 
O=AUTO STOICH. RATIO OF CARBON OXIDES 1=USER ENTRY OF RATIO 
BOUNDARY CONDITION PARAMETER 

1=AUTO LINEAR TEMPG, TEMPA, CONC. O02 O=USER ENTRY (SEE BELOW) 
LOCATION OF HEATER O=HEATER AT ZERO 1=HEATER AT XLENGTH 


500 NSTOP MAX NUMBER OF TIMES ALLOWED THROUGH INTEGRATION LOOPS 
0 ITRIP LEAVE ZERO UNLESS RESTARTING PROGRAM IN SURF RECESSION MODE 


--GRAPHITE TEMP NODE 1]*THESE 3 ITEMS ARE REQUIRED ONLY 
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800.E0 


.0172E0 


0.DO 
.001D0 
QO 5D0 
-O1ED0 
2025D0 
<05D0 
. 1D0 
“LS B0 
.2D0 
22500 
73D0 
sooU0 
.4D0 
.45D0 
1550 
,5 500 
.6D0 
7c 5 D0 
27 D0 
.738D0 
,3D0 
soo D0 
.9D0 
79500 


- AIR TEMP NODE, 1---] ]*IF IIC=1. THEY MUST FOLLOW "ITRie eee 
--OXYGEN CONC NODE 1--]*MANUAL ENTRY OF NODE LOCATIONS IF USED 
*xex*x**THE FOLLOWING DATA IS THE X LOCATION OF THE NODES 
*THESE LOCATIONS ARE ONLY NEEDED IF IELEM=1 
*REQUIRING NODAL LOCATION INPUT. IF IELEM=0 
*THE COMPUTER CALCULATES A LINEAR SPACING OF THE 


*NODES. 

* 

* 

* IF THIS MANUAL ENTRY IS USED IT MUST 
* FOLLOW INPUT OF THE VARIABLE "ITRIP" 
* AND HAVE A VALUE FOR EACH NODAL POINT 
* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


* 


1.DO *******END OF NODAL POINT LOCATION INPUT * * ® # i ti i i ee ei 


800.D0 
800 .D0 
800.D0 
800.D0 
800.D0 
800 .D0 
800 .D0 
800.D0 
800.D0 
800.D0 
800 .D0 
800.D0 
800 .D0 


SSSSSSSSSTART OF MANUAL ENTRY OF GRAPHITE TEMP FOR EACH NODE 


SSSS$THIS INFORMATION ON TEMPERATURE AND CONCENTRATION 
IS ONLY REQUIRED IF IIC=0 REQUIRING MANUAL DATA INPUT 
FOR EACH NODAL POINT. THE $SS$$S$$ DATA MUST FOLLOW DATA 
FOR "ITRIP" IF NODAL LOCATIONS ARE NOT BEING MANUALLY 
INPUT (IF IELEM=1)OR ELSE THIS DATA FOLLOWS NODAL POINT 
MANUAL ENTRY LOCATIONS 


DATA MUST BE GIVEN FOR EACH NODAL POINT 


TLAATATATAMMMNAUMUMUNAMNUMUUAMUMNUNOUMAUNUNUNURAMNN 


SSSSSSLAST MANUAL ENTRY OF GRAPHITE TEMPERATURES%S%S%S%% 
SSSSSSSTART MANUAL ENTRY OF AIR TEMPERATURE EACH NODE## 


tt Tr 


NNN Ua 
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ATAU MUNNONUNUNUNUNYN 


SSSSLAST MANUAL ENTRY OF AIR TEMPERATURE EACH NODE 
SSSSSTART MANUAL ENTRY OF OXYGEN CONCENTRATION EACH NODE 


$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 


SSSSLAST MANUAL ENTRY OF OXYGEN CONCENTRATION EACH NODE 
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APPENDIX E: 


31 NODAL POINTS 
AT TIME= 
O0018S 


MEDIUM THICKNESS= 


.4166666666667D-01 


NP 


OO O°O Oo: © Oo OO © oc Co Ooo OO © oo @ @ oO © © 0 0 02 ea 6 


X 


.O0O00000000000000E+00 
.1388888888888889E-02 
wZITIVIIIT IF ITT 77 eeeO2 
-4166666666666666E-02 
5555555555555 5 55802 
.6944444444444443E-02 
pO3333953233253 318-02 
-97Zecceeececcecin- OZ 
op er Wes ea Wes Wo a a op a Ba Ee Ro) ye oa 
.1250000000000000E-01 
-1388888888888889E-01 
S15 277977720777 77 -On 
-L6E66666666666666E-01 
~£8055555555555552-02 
.1944444444444444E-Q01 
~2083335333333333E-01 
pee2e2Z2222222222228-01 
2361417 Didier EO 
RA cas ys Woke belles) o se iss spars opk 
.2638888888888888E-01 
SQITIQVIPICII 197 ae Ok 
-2916666666666666E-01 
JOSS 5555555 555 5.55)- Orb 
.3194444444444444E-01 
733393555399 393538- 08 
»9472222222222222E-01 
p36 TPIT IAAT IIa OE - OL 
-o 749999999999 999 8-0. 
-3888888888888888E-01 
4027777 II i777 OL 
-4166666666666666E-01 


oOo © 0 © OO © oo CO ©] © OC CO eG Oo © Oo 6 © 6 6 © 6 6 2 & 


.4166666666667D-01 


TEMP -G 


-4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
.4000000000000000E+03 
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SAMPLE OUTPUT DATA FROM COMBUSTION CODE 


DIAMETER 
.6944444444444444E-03 
.6944444444444444FE-03 
.6944444444444444E-03 
.6944444444444444E-03 
-6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444FE-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
-6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 
-6944444444444444E-03 
.6944444444444444E-03 
.6944444444444444E-03 


AT TIME= 


SoS 


MEDIUM THICKNESS= 


.4166666666667D-01 


NP 


NO F 


OO OOO 0 © CC Oo © OC @ OCG OO OO eG @ OC OG 6 40 2 oo 


x 


.0O000000000000000E+00 
-1388888888888889E-02 
P2799 777777777177 718E-02 
.4166666666666666E-02 
S555 5555555555595 - 02 
.6944444444444443E-02 
meses sseseoS5oSS531E-02Z 
Moecezeeeaa2e2c2LE- 02 
eed eee lib On 
-1250000000000000E-01 
-1388888888888889E-01 
USP 9 7) 1607 Tat A A a ag ad fe Ol 
-1666666666666666E-01 
Pio Sp 55555555 555E 08 
.1944444444444444E-01 
Ne0GS2555555333332E-01 
weaeegageaaceecaean-OL 
resol et tle Pe - Or 
pea 9999999999999 E-o0L 
-2638888888888888E-01 
UZIP IIo aI TAIT ITE O01 
-2916666666666666E-01 
 USs5 505555555 55E-02 
.3194444444444444E-01 
SooeoSo 2665 5553333E-01 
Sag ae cee Lala OL 
Soler ttt OR -02 
Bo ages IO IIo II999E-O1l 
-3888888888888888E-01 
§rlUZI77 777 Pig TIT TIE =0L 
.4166666666666666E-01 


Oo ooqo0oqcoo coocnoocooodlomdmhmd lOO OO ee ea eae © eo eae 


.4166666666667D-01 


TEMP-G 


-6171655414743030E+03 
.5781984237740606E+03 
*55 121344 753523755E+03 
.5261308274882458E+03 
POs 56fOR36879362E+03 
.4836232366654954E+03 
StGosoo5 79 55059Z35E403 
.4517249092494148E+03 
-4395833460910138E+03 
-4297262006976060E+03 
.4218950356130471E+03 
.4158052068913094E+03 
.4111679694792949E+03 
.4077084937514894E+03 
-4051785413371006E+03 
-4033636888937181E+03 
.4020858170994085E+03 
PtOtZOzOorT I Sls2163sE+03 
.4006011841982683E+03 
.4001993564763142E+03 
.3999347556308410E+03 
.3997630187488629E+03 
-3996530135102078E+03 
13995833611737365E+03 
~Ss9a 552 be22550029E+03 
-3995124994579232E+03 
.3994956957918823E+03 
-3994854173856393E+03 
.3994793187448693E+03 
.3994760668482828E+03 
-3994750372395511E+03 
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oooooooocoeoooo°o°o°ooe°oeoeo0°o0 0 0coe0°0 0 Co 0co0co0Cco°c °o 


DIAMETER 
.6944444021540727E-03 
-6944444306055943E-03 
.6944444380346299E-03 
-.6944444413489135E-03 
-6944444428425697E-03 
-6944444435425437E-03 
-6944444438859773E-03 
-6944444440630632E-03 
-6944444441590080E-03 
-6944444442134259E-03 
-6944444442455372E-03 
-6944444442651040E-03 
-6944444442773239E-03 
-6944444442850919E-03 
-6944444442900914E-03 
-6944444442933369E-03 
.6944444442954588E-03 
.6944444442968574E-03 
.6944444442977906E-03 
-6944444442984252E-03 
-6944444442988702E-03 
.6944444442991956E-03 
.6944444442994469E-03 
-6944444442996526E-03 
-6944444442998301E-03 
-6944444442999913E-03 
-6944444443001414E-03 
.6944444443002882E-03 
.6944444443004240E-03 
.6944444443005773E-03 
.6944444443006692E-03 


AT TIME= 


dW fare = a eer 


MEDIUM THICKNESS= 


.4166666666667D-01 


NP 


NO 


Oo oO COCO OO OCOOCOCOO COO OOOO OO Oo eo eo eo oa & 


x 


.0000000000000000E+00 
-1388888888888889E-02 
S207 ITI Id 171297 1 Be -O7 
-4166666666666666E- 02 
[5555 5555550555555 a0 of 
.6944444444444443E-02 
PO 3335 3955533333 1R 02 
AS Ps PAPA PAP APAPA Pt OOP A PAPAIN Th 0) 2 
ips ia 2 a Vga bg Wega LL Ba Wp i a Up (9 
-L250000000000000E-01 
.1388888888888889E-01 
SES 2e 77 V7 oe 
.L666666666666666E-01 
PI60SS555555555554- 08 
.1944444444444444E-01 
2209035533555533335-0 
PAPI DLP PAPAP PAP IEDAPAED ISI AND SI A 
£235 607 Tae Oe 
:2499999999999999E—02 
-.2638888888888888E-01 
ATI TTI tee ee Od 
.2916666666666666E-01 
-3 055555555555 5555-01 
.3194444444444444E-Q1 
39955333 35353533335»-01 
se4 7 2222222222222 -01 
ee syialab ala babs bal alaaba ga 02; <0 al 
33 1499999999999 99> 01 
.3888888888888888E-01 
4027777777 771777E-01 
-4166666666666666E-01 


0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 


.4166666666667D-01 


TEMP-G 


.6944483272376198E+03 
.6582952483489812E+03 
.6330284155199 859E+03 
.6085590421656422E+03 
.5853068872440880E+03 
.56339151Z227869359E+08 
.5429230303669187E+03 
.5239878470762538E+03 
.5066449982116862E+03 
.4909237200244683E+03 
.4768226385980613E+03 
.4643106220267744E+03 
.4533292103391948E+03 
.44379638472023 81E+03 
.4356113128753018E+03 
.4286596365129955E+03 
.4228188585575783E+03 
.4179634359788205E+03 
.4139692729896838E+03 
.4107174182763530E+03 
.4080968793508916E+03 
-.4060065623126001E+03 
.4043564178987307E+03 
.4030679223461711E+03 
.4020740463757375E+03 
.4013188723800529E+03 
.4007570137582123E+03 
.4003529789248906E+03 
.4000805996527675E+03 
- 3999226594325) 25E+03 
-3998706998552636E+03 
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DIAMETER 
.6944440200051797E-03 
.6944442750792861E-03 
-6944443547335164E-03 
.6944443969186038E-03 
.6944444188805062E-03 
.6944444303859835E-03 
.6944444364794849E-03 
.6944444397590113E-03 
.6944444415612797E-03 
.6944444425767153E-03 
.6944444431650298E-03 
-6944444435161306E-03 
-6944444437320398E-03 
-6944444438687254E-03 
-6944444439576309E-03 
-.6944444440168827E-03 
-6944444440572177E-03 
-6944444440851723E-03 
.6944444441048360E-03 
.6944444441188342E-03 
.6944444441288955E-03 
.6944444441361831E-03 
.6944444441414949E-03 
-6944444441453872E-03 
.6944444441482521E-03 
.6944444441503701E-03 
.6944444441519371E-03 
-6944444441530992E-03 
-6944444441539341E-03 
.6944444441545535E-03 
.6944444441548335E-03 


AT TIME= 


262 565 


MEDIUM THICKNESS= 


.4166666666667D-01 


WOODMTNHNU PP WHF 
UO 


oo oO 0 0 OOo 0 OO SO © OO © O20 Oe 2 2 6 ea oo 


Xx 


.O0O00000000000000E+00 
.1388888888888889E-02 
(2) Te 1 te a ie Oe — OZ 
-4166666666666666E-02 
poo SS Sooo Soo SoR-O2 
.6944444444444443E-02 
fees 5555565333318 -02 
5) IE PAPA PAOD PLEA PLIL PAPA Y Das OW 
rll Ee TT: TET TT Fe a Bs Bo a Bp Ona 8 
.L250000000000000E-01 
.1388888888888889E-01 
SES 2) i 7 7 TE OL 
. 1L6E66666666666666E-01 
iO > Seeias 555555 58-> 01 
.1944444444444444E-01 
P20833333333333533E-01 
p2222e22er~-2ce222c48-O01 
resol vite-O1 
12497999 999999999E-01 
.2638888888888888E-01 
21 egg T7771 7TE- OL 
.2916666666666666E-01 
(JO SSSSS555555555E-01 
.3194444444444444E-01 
po oS oo Sooo SSSSE-01 
~3472222222222222E-01 
oo cleaner Pie 0 1 
Poe got III 99999 99E-01 
.3888888888888888E-01 
P2O2Z777 77771717777 1E-O1 
-4166666666666666E-01 


0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 


.4166666666667D-01 


TEMP-G 


.6306595038438272E+03 
.6399241320641830E+03 
.6447192511167083E+03 
-6406985836760575E+03 
.6302556153620338E+03 
.6156188718308147E+03 
.5986698701861527E+03 
.5807979040862021E+03 
.5629290453016817E+03 
Poaesoe25613672719E+03 
moog 2950993824 7957E+03 
.5139161110240054E+03 
.4997594668844450E+03 
-4867985002414993E+03 
.4750316490868472E+03 
.4644343493352758E+03 
-4549653630104841E+03 
-4465709326513697E+03 
-4391880102692778E+03 
-4327471084903999E+03 
.4271749671005763E+03 
42239 70689333552E+03 
-4183399816653412E+03 
.4149334943623446E+03 
a bos '5457292E+03 
.4098188488450426E+03 
-4080025338059548E+03 
.4066233572360078E+03 
.4056520065572266E+03 
-4050713541567405E+03 
.4048773671370597E+03 
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DIAMETER 
.6944429267898032E-03 
.6944437360540249E-03 
.6944440181040645E-03 
.6944441916810881E-03 
.6944442957791966E-03 
.6944443573613218E-03 
.6944443933015021E-03 
.6944444141127005E-03 
-6944444261547878E-03 
.6944444331631678E-03 
.6944444372871645E-03 
.6944444397507130E-03 
.6944444412492834E-03 
.6944444421795931E-03 
.6944444427698287E-03 
.6944444431527831E-03 
.6944444434068504E-03 
.6944444435790780E-03 
.6944444436982122E-03 
-6944444437821565E-03 
.6944444438422836E-03 
.6944444438859626E-03 
.6944444439180605E-03 
.6944444439418518E-03 
-6944444439595736E-03 
.6944444439727761E-03 
.6944444439825334E-03 
.6944444439895999E-03 
.6944444439944511E-03 
.6944444439974719E-03 
.6944444439985787E-03 
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